#include "n2.h"
#include <iostream>
#include <string>
#include <sstream>
#include <cmath>
#include <ctime>
#include <cstdio>
#include <cfloat>
#include <boost/program_options.hpp>



void printVec(float *data, unsigned texSize)
{
	unsigned nData = texSize * texSize * 4;

	for (unsigned i = 0; i < nData; i += 4)
		std::cout 
			<< "(" << data[i] << ","
			<< data[i+1] << ","
			<< data[i+2] << ","
			<< data[i+3] << ")"
			<< std::endl;
}
void plasterVec(float *data, unsigned texSize, 
				float region, float mass)
{
	unsigned nData = texSize * texSize * 4;
	srand(0);
	for (unsigned i = 0; i < nData; i+=4) {
		//data[i] = (float) rand() / RAND_MAX;
		data[i] = (double) (rand() - RAND_MAX / 2) 
			/ RAND_MAX * region;
		data[i+1] = (double) (rand() - RAND_MAX / 2) 
			/ RAND_MAX * region;
		data[i+2] = (double) (rand() - RAND_MAX / 2) 
			/ RAND_MAX * region;
		data[i+3] = mass;
	}
}
void describeError(float *fa, float *fb, unsigned texSize)
{
	float d, dMax = 0;
	double dSum = 0;
	unsigned iMax = 0;
	unsigned nData = texSize * texSize * 4;
	for (unsigned i = 0; i < nData; i+=4) {
		if ((d = abs(fa[i] - fb[i])) > dMax) {
			dMax = d;
			iMax = i;
		}
		dSum += d;
		if ((d = abs(fa[i+1] - fb[i+1])) > dMax) {
			dMax = d;
			iMax = i;
		}
		dSum += d;
		if ((d = abs(fa[i+2] - fb[i+2])) > dMax) {
			dMax = d;
			iMax = i;
		}
		dSum += d;
	}
	std::cout
		<< "Error: max=" << dMax 
		<< ", avg=" << dSum / nData
		<< ", (" << fa[iMax] << " " << fa[iMax+1] << " " << fa[iMax+2]
		<< ") (" << fb[iMax] << " " << fb[iMax+1] << " " << fb[iMax+2]
		<< ") at " << iMax
		<< std::endl;
}
void computeOnCPU(float *bodies, float *forces, 
				  float eps, unsigned texSize)
{
	double r2, r3, rij[3];
	unsigned nData = texSize * texSize * 4;
	double tmp[3];


	for (unsigned i = 0; i < nData; i+=4) {
		//forces[i] = forces[i+1] = forces[i+2] = forces[i+3] = 0;
		tmp[0] = tmp[1] = tmp[2] = 0.0f;
		for (unsigned j = 0; j < nData; j+=4) {
			rij[0] = bodies[j] - bodies[i];
			rij[1] = bodies[j+1] - bodies[i+1];
			rij[2] = bodies[j+2] - bodies[i+2];
			r2 =  rij[0] * rij[0]
				+ rij[1] * rij[1]
				+ rij[2] * rij[2]
				+ eps;
			r3 = r2 * sqrt(r2);
			/*
			forces[i] += bodies[j+3] * rij[0] / r3;
			forces[i+1] += bodies[j+3] * rij[1] / r3;
			forces[i+2] += bodies[j+3] * rij[2] / r3;
			*/
			tmp[0] += bodies[j+3] * rij[0] / r3;
			tmp[1] += bodies[j+3] * rij[1] / r3;
			tmp[2] += bodies[j+3] * rij[2] / r3;
		}
		/*
		forces[i] *= bodies[i+3];
		forces[i+1] *= bodies[i+3];
		forces[i+2] *= bodies[i+3];
		*/
		forces[i] = (float) (tmp[0] * bodies[i+3]);
		forces[i+1] = (float) (tmp[1] * bodies[i+3]);
		forces[i+2] = (float) (tmp[2] * bodies[i+3]);
	}
}
namespace po = boost::program_options;
int main(int argc, char **argv)
{
	time_t start, stop;
	unsigned texSize;
	float eps, mass, region;
	std::string sourceName;

	po::options_description opt("Options");
	opt.add_options()
		("show_help",
		"Show this help message")
		("tex_size",
		po::value<unsigned>(),
		"The square size")
		("dist_eps",
		po::value<float>()->default_value(1.0e-6f),
		"Minimum distance of bodies")
		("mass",
		po::value<float>()->default_value(10000.0f),
		"Default body mass")
		("region",
		po::value<float>()->default_value(10000.0f),
		"Default body region")
		("cg_source",
		po::value<std::string>()->default_value("shader.cg"),
		"Shader source file name");
	po::variables_map vm;
	po::store(parse_command_line(argc, argv, opt), vm);
	po::notify(vm);

	if (vm.count("show_help") || !vm.count("tex_size")) {
		std::cout << opt << std::endl;
		return 0;
	}
	texSize = vm["tex_size"].as<unsigned>();
	eps = vm["dist_eps"].as<float>();
	sourceName = vm["cg_source"].as<std::string>();
	mass = vm["mass"].as<float>();
	region = vm["region"].as<float>();
	
	unsigned nData = texSize * texSize * 4;
	std::cout
		<< "texSize="
		<< texSize
		<< ", nData="
		<< nData
		<< std::endl;

	float *bodies = new float[nData];
	float *fgpu = new float[nData];
	float *fcpu = new float[nData];
	plasterVec(bodies, texSize, region, mass);
	//printVec(bodies, texSize);

	GPGPU gpgpu(texSize);
	gpgpu.init(argc, argv);
	Kernel kern(texSize, sourceName, rect_nv_rgba_32);
	kern.init();

	//for (unsigned i = 0; i < 100; i++)
	time(&start);
	kern.calc(bodies, fgpu, eps);
	time(&stop);
	std::cout
		<< "Time took on GPU [sec]: "
		<< difftime(stop, start)
		<< std::endl;
	//printVec(fgpu, texSize);

	time(&start);
	computeOnCPU(bodies, fcpu, eps, texSize);
	time(&stop);
	std::cout
		<< "Time took on CPU [sec]: "
		<< difftime(stop, start)
		<< std::endl;
	//printVec(fcpu, texSize);

	describeError(fgpu, fcpu, texSize);

	return 0;
}
